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Abstract 

By using the properties of orthogonal polynomials, we present an exact 
unitary transformation that maps the Hamiltonian of a quantum system 
coupled linearly to a continuum of bosonic or fermionic modes to a Hamil- 
tonian that describes a one-dimensional chain with only nearest-neighbour 
interactions. This analytical transformation predicts a simple set of rela- 
tions between the parameters of the chain and the recurrence coefficients 
of the orthogonal polynomials used in the transformation, and allows the 
chain parameters to be computed using numerically stable algorithms that 
have been developed to compute recurrence coefficients. We then prove 
some general properties of this chain system for a wide range of spectral 
functions, and give examples drawn from physical systems where exact 
analytic expressions for the chain properties can be obtained. Crucially, 
the short range interactions of the effective chain system permits these 
open quantum systems to be efficiently simulated by the density matrix 
renormalization group methods. 



1 Introduction 

All quantum systems encountered in nature experience random perturbations 
due to their coupling to degrees of freedom of their local environments. Informa- 
tion about these degrees of freedom are not normally accessible, and to correctly 
predict the results of experiments, these degrees of freedom must be averaged 
over. This averaging introduces qualitatively new features into the otherwise 
unitary dynamics of the quantum system, and typically induces an effectively 
irreversible dynamics which drives the system towards an equilibrium with its 
environment [I]. For quantum systems, this evolution towards equilibrium not 
only involves energy transfer to the environment, but can also cause the loss of 
coherence in the system state, often on a much faster timescale than the energy 
relaxation [3 [5] . This latter process of decoherence rapidly destroys quantum 
mechanical effects arising from the existence of phase coherence in the state of 
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the system, and must always be considered when analyzing the results of ex- 
perimental investigations of quantum phenomena. Decoherence also presents a 
major problem for the emerging field of quantum information and technology, 
a branch of physics that aims to directly harness quantum mechanical effects to 
create novel and highly effective devices which could greatly outperform their 
classical counterparts. This latter issue has attracted considerable interest in 
recent years, and the study of quantum systems in contact with environments, 
often referred to as 'open-quantum systems', is a key problem in contemporary 
physics [2]. 

The theoretical description of a quantum system coupled to environmental 
degrees of freedom has been developed over many decades, and of particular 
important in these studies is the spin-boson model (SBM) [U [3], one of the 
simplest non-trivial models of open-system dynamics. This model describes a 
single two-level system (TLS) coupled linearly to the coordinates of an environ- 
ment consisting of a continuum of harmonic oscillators. Despite its simplicity, 
this model cannot be solved exactly and shows a rich array of non-Markovian 
dynamical phenomena, including a type of quantum phase transition between 
dynamically localized and delocalized states of the TLS for so-called Ohmic and 
sub-Ohmic baths at zero temperature [21 HI HI IB ISl H] ■ These interesting dy- 
namics normally appear in situations where the effective interaction between the 
TLS and the oscillators cannot be treated perturbatively, and although many 
analytical and formal approaches have shed considerable light on the nature of 
the effects, the explicit time evolution of the TLS in these regimes must nor- 
mally be simulated numerically. These simulations can be loosely divided into 
methods in which the environmental degrees of freedom are eliminated in order 
to derive an effective equation of motion for the TLS, and those in which the 
complete many-body dynamics of the TLS and the environment are computed. 
The former approach includes the hierarchy approach and various numerical 
path integral techniques, while the latter includes exact diagonalization and the 
numerical reorganization group (NRG) approaches jTH! . 

Recently, the time-adaptive density matrix renormalisation group (t-DMRG) 
[Tl] has also been added to this category. This approach provides numerically- 
exact simulations for the dynamics of one-dimensional systems with short-range 
interactions, and is regarded as one of the most powerful and versatile methods 
in condensed matter, atomic, and optical physics for the study of strongly- 
correlated many-body quantum states. The recent study of Prior et al. [S] 
successfully used a t-DMRG approach for a two-spin system where each spin 
was in contact with an oscillator environment with arbitrarily complex spectral 
densities. The implementation and accuracy of this approach rests on an exact 
mapping between the canonical SBM and a one-dimensional harmonic chain 
with short range interactions. The reasons why this form of Hamiltonian allows 
the t-DMRG to efficiently simulate extremely large chains is disscussed at length 
in Ref. TT . 

The idea of this mapping was first introduced in the NRG studies of Bulla 
et al. [9 , who first discretise the continuous environment and then perform 
a recursive numerical technique to bring the truncated finite Hamiltonian to 
the desired chain form. The dynamics of this chain can then be numiercally 
simulated [12]. In this paper we show how the formal properties of orthogonal 
polynomials can be used to implement an exact unitary transformation which 
maps an arbitrary continuous SBM directly onto the nearest-neighbour 1 — 
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D chain without any need for discretisation. In several important cases we 
can actuahy perform the mapping analytically and predict the energies and 
couplings of the chain in closed form. Our analyisis of the mapping shows that 
these energies and couplings are simply related to the recurrence coefficients 
of orthogonal polynomials, and when analytical results are not available, we 
can make use of sophisticated numerical techniques that compute recurrence 
coefficients with high precision and which do not suffer from the numerical 
instabilities which often appear in the standard recursive technique. Not only 
does this transformation aid the simulation of non-perturbative open-system 
dynamics, we shall show that the use of orthogonal polynomials also reveals 
some universal features of TLS dynamics which are independent of the spectral 
density of the environmental oscillators, and which may point the way to even 
more efficient simulations of strongly non-markovian open-system dynamics. 

The plan of this paper is as follows. In section [5] we introduce the formal 
properties of orthogonal polynomials needed to implement the chain mapping, 
and in section [3] we introduce the general model of an open-quantum system 
which can be mapped onto &. 1 — D chain suitable for t-DMRG simulation. 
This section also contains details about the transformation itself and analyt- 
ical expressions for the frequencies and nearest-neighbour interactions of the 
chain are presented and discussed. Section 2] gives explicit examples of this 
transformation for the SBM, and 14.21 deals with an extension of the orthogonal 
polynomial method to systems interacting with a discrete number of environ- 
mental modes. These discrete environments may be found in certain physical 
systems, and are also important in the NRG description of open-quantum sys- 
tems where the continuous spectrum of the environment is approximated by a 
discrete set of oscillators [9] . We then end with some general conclusions and a 
brief recapitulation of the main results. 

2 Orthogonal Polynomials 

In this section we present a few useful properties of orthogonal polynomials that 
will be employed in the following. This section also serves to fix our notation. 

2.1 Basic concepts 

Let P be the space of the polynomials with real coefficients (the following may 
also be extended to complex, but real coefficients will be enough for our pro- 
poses), and P„ C P the space of polynomials of degree equal to n, this is 



where the G R and a; is a real variable. 

Definition 2.1. A polynomial Pn{x) G Pn is called monic if its leading co- 
efficient satisfies a„ ~ 1, will be denoted by Pnix) = 7r„(a;). Of course every 
polynomial of P„ can become monic just dividing it by a„. 

Definition 2.2. Let ^{x) be a nondecreasing and differentiable function on 
some interval [a, b] € R, and assume that the induced positive measure dfiix) 



Pn{x) £ P„: 



n 




m=0 
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has finite moments of all orders: 

l-b 



/ 

J a 



x^diJ,{x) < oo, r = 0, 1, 2, . . . 
Then for any p{x), q{x) e P we define an inner product as 

{P,q)ij. = / p{x)q{x)dii{x). 

J a 



Of course it defines a norm via 

\ 1/2 

p'^{x)dii{x) ) > 0, 



and satisfies the Cauchy-Schwarz's inequality, 

\{p,q),\<\\pUq\\,. 

Definition 2.3. Two polynomials p{x), q{x) G P are said to be orthogonal with 
respect to the measure d^{x) if its inner product is zero (p, q)^^ = 0. 

Definition 2.4. A set of monic polynomials {7r„(x) € Pn,?^ = 0,1,2,...} is 
called an orthogonal set with respect to the measure dfi{x) if 

(TTn, Trm)fi = ||7r„||^^„,m for n, TO = 0, 1, 2, . . . 

Sometimes it is useful to normalize the monic polynomials, Pn{x) = 7r„(a;)/|l7r„(.x)||^, 
such that {pniPm)n = <5ra,m) then the set is called orthonormal, however note 
that in general the normalized polynomials lose the monic character. 

Theorem 2.1. For any (non-zero) measure dfj.{x) there exist a unique family 
of monic orthogonal polynomials {7r„(a;) e Pn, n = 0, 1, 2, . . .}. 

Proof. The proof is by construction. Applying the well-known Gram-Schmidt 
orthogonalization procedure to the sequence {1, x, x^, . . .}, the monic orthogonal 
polynomials are constructed from ttq = 1 using the formula: 

fe-i 



TTkix) = mk{x) - > ^ TTn 



with mk{x) = x'^. □ 

Theorem 2.2. The set of monic polynomials {7r„(x) G P„, n = 0, 1, 2, . . .} is a 
basis of the space P. More specifically if p{x) is a polynomial of degree less or 
equal to r, then it can be uniquely represented by 

r 

P{x) = '^Cn7rn{x) (1) 

for some real constants Ck ■ 
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Proof. The proof is by mathematical induction. If r = the proof is trivial. 
Assume that the statement is true for r — 1, then there is a unique Cr such 
that p(x) — CrTTrix) = Pr-i{x) is polynomial with degree less or equal to r — 1. 
Moreover, if the set of monic polynomials is orthogonal, taking the inner product 
on ID) with respect to Trm{x) we obtain 




There are of course many more properties like these which are similar to 
those found in linear algebra spaces; however we are not going to go further 
in this way here, though the reader can find more information in references 
[131 1141 115) for example. From now on we omit the subindex ^ in the inner 
product and norms in order to make the notation less dense. 



2.2 Recurrence Relations 

One of the most useful properties for our interests is the recurrence formula of 
the orthogonal polynomials. Here we divide the section into two cases, monic 
and orthonormal polynomials. 

Theorem 2.3. Let {7r„(a;) G P„, n = 0, 1, 2, . . .} be the set of monic orthogonal 
polynomials with respect to some measure dfi{x), then 

TTk+iix) = {x - ak)Trk{x) - iSkTTk-iix), fc = 0,1,2... (2) 

where t:^i(x) = 0, and the recurrence coefficients are: 

ttfc = -, ^, Pk = 7 r- (o) 

Proof. By theorem [2.21 Tr^+i (x) — xiTk{x) is a polynomial of degree less or equal 
to /c, we can write it in terms of monic polynomials of degree less or equal k: 

k 

■Kk+i{x) - xTTkix) = ^ c„7r„(a;), 
and since the set is orthogonal 

(TTfc+l - CTfe, TTfc) _ {xTTk,TTk) 



Ck 



similarly 

(tt/c+i - XTTfc, 7rfc_i) _ {xnk,Trk-i) _ {T^k, xiTk-i) _ (7rfe,7rfc) 



(TTfe-i, TTfe-i) (7rfc_i, 7rfe_i) (tt/c-i, 7rfe_i) (7rfe_i, 7rfc_i) 



where in the last step we have used the fact that terms of degree less than k 
are orthogonal to 7rfe(x), because they can be expressed as linear combination 
of {nn{x) £ P„, n = 0, 1, 2, . . . , fc — 1} which are all orthogonal to Trk{x). By the 
same argument, it is easy to see that the remaining coefficients are zero c„ = 
for n = 0, . . . , fc — 2. So by comparing with equation ^ we obtain ak — —Ck 
and = -Cfc-i. □ 
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Theorem 2.4. Assume that we have a set of orthonormal polynomials {pn{x) £ 
P„, n = 0, 1, 2, . . .} then, the recurrence relation is: 

Pk+i{x) ^ {CkX - Ak)pk{x) - BkPk-i{x), fc = 0,1,2... (4) 

where p^i{x) = 0. If we denote by k„ the leading coefficient of pn{x) and by k!^ 
the second leading coefficient, i.e. Pn{x) = Knx"^ + k'^x"^^ + . . we get 

Ak = ^(^-^), Bk^^^^±^, = ^ (5) 



Proof. The proof is not more complicated than in the monic case, if we chose 
to /c, in the same fashion the choice 



Cfe = then, pk+i{x) — Ckxpk{x) is a polynomial of degree lower of equal 



Ckn'k ''^'k+l l^k+l ( Kfc '*fe+l 
Ak = = 



Kk Kk \Kk Hk+1 

makes pk+i{x) — (CkX — Ak)pkix) a polynomial of degree lower of equal to k—1. 
By orthogonality the inner product of it with p„ for n < k — 2 are zero, and the 
inner product with pk-i yields 

{pk-i,Pk+i - {Ckxpk + Ak)pk) = -Ck{pk-i,xpk) = -Ck{xpk~i,Pk) 

^ Hk-l i~ ~ > ^ Hk^l 

= -Ck {Pk,Pk) = -Ck , 

Kfc Kk 

SO by taking into account that this is the coefficient which goes with Pk-i, 
replacing the value of Ck gives the value of in ([5|) . □ 

The relation between both cases is given by the next proposition. 

Proposition 2.1. Let {7r„(a;) e P„, n = 0, 1, 2, . . .} a set of monic orthogonal 
polynomials and pn{x) = 7r„(a;)/||7r„(x)|| their orthonomalized version, then the 
coefficients for the respective recurrence relations are related by: 



Ak^^^, Bk^J^ andCk = ^=. (6) 

Proof. Inserting 7r„(a;) = ||7r„(x)||p„(x) in ^ and dividing by |l7rfc-|_i|| one ob- 
tains 

Pk^\\X) = -u l—^M ~ ^kjPk{x) - J: —--Pk-l{x). 

Employing f3k — ||7rfc|p/||7rfc_i||^ and comparing with the above (jj]) we obtain 

(El). □ 

2.3 Boundness properties of the recurrence coefficients 

A particular important property for our work will be the behaviour of the monic 
recurrence coefficients a„ and /3„ when n is increasing. To start, we have the 
following theorem. 
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Theorem 2.5. Let the support interval [a, b] of dfj,{x) be finite, then for each 
k = 0,1,2,... 

a < ak <b, 

Proof. The first relation follows immediately from ak = ^^^^'^^-j because tt^ (cc) 
is positive for every x and a < x < b inside of the integration domain. For the 
second relation we have by orthogonality 

(7rfc,7rfe) = (x7rfc_i,7rfc), 

and the Cauchy-Schwarz's inequality gives 

(TTk^T^k) < |a;|||7rfc_i||||7rfc|| < max{|a|, |6|}||7rfc_i||||7rfc||, 

finally the division of both sides by ||7rfe|| and ||7rfc„i|| and squaring leads straight- 
forwardly to the second bound. □ 

Note that for measures with unbounded support, the monic recurrence co- 
efficients are not bounded as n increases; see for example the case of Laguerre 
and Hermite polynomials US'. However, for measures with bounded support 
further useful properties about the recurrence coefficients can be obtained, to 
do this we need the following results. 

Proposition 2.2. Let {7r„(a;) G P,i, n = 0, 1, 2, . . .} be a set of monic orthonor- 
mal polynomials with respect to some measure d^{x) supported in [—1, 1], with 
recursion coefficients ak and /3k defined according to (0)- If we change the sup- 
port of the measure to [a, b] by means of a linear transformation y ~ mx + c, 
where 

b — a a + b , . 

"^=^, c=^, (7) 

the recursion coefficients a'k and of the transformed monic polynomials {7r^(y) 
G P„, n = 0, 1, 2, . . .} are given by 

a'k = mak + c, jS'k ^ m^jSk. (8) 

Proof. The map y — mx + c transform the measure dfi{y) ~ md[i(x) and every 
monic polynomial i^kix) — 7rfc(^^) which is not a monic polynomial in To 
make these polynomials monic in we need to multiply by , that is, the set 
of transformed monic polynomials is {7r^(?/) = m^Trk{x) G P„, n = 0, 1, 2, . . .}. 
Then we apply the definitions Q 

Ct 7,, 



j^2k+l {jni^x-Kk , TTfe ) -I- c{-Kk , TTfc ) ) 



and similarly 



i2fe+i(7rfe,7rfe) 



T'k^TT'k) _ TO2^-+l(^fc,7rfe) _ 2 



mak + c, 



K-i^k-i) m2(fc-i)+i(7r,_i,^,..i) 



m'Pk- 

□ 
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Definition 2.5. A measure dij.{x) = w{x)dx supported in [—1,1] belongs to 
the Szego class if 

lnu;(a;) , 



-dx > —oo. 



Theorem 2.6 (Szego). Let {p„(x) £ Pn,n ~ 0,1,2,...} orthonormal polyno- 
mials with respect to some Szego class measure dfi{x) = w{x)dx, then in the 
asymptotic limit n ^ oo we have the approximate expressions for the first lead- 
ing coefficient k„ and for the second one k^, 



: exp 



r3/2 



1 lnw{x) 
^ xlnwlx) 



dx 



1 vr^~^2 



rdx 



exp 



1 

2^ 



^ lnw{x) 

1 



dx 



(10) 

(11) 
□ 



Proof. The proof can be founded in chapter XII of 15 . 
Proposition 2.3. The coefficients An, Bn and Cn have the asymptotic values 
lim A„ = 0, hni B„ = 1, lim C„ = 2. (12) 

Proof. It follows immediately from ^ and the previous result ((TO)) and dTTI) . □ 

The speed of convergence of An, Bn and C„ depends on the measure in 
a non-trivial fashion (see details [H]). We will demonstrate this explicitly in 
equations (P7| and (pS)) for specific spectral densities. 

Corollary 2.1. Let {7r„(x) G P„, n = 0, 1, 2, . . .} be a set of orthogonal monic 
polynomials with respect to some measure d^{x) which belongs to the Szego 
class, then their recurrence coefficients have the asymptotic values 



lim q;„ = 0, lim /3„ 



1 

4' 



Proof. See proposition 



□ 



Proposition 2.4. Let {7r„(a;) e F„, n = 0, 1, 2, . . .} he a set oforthogonal monic 
polynomials with respect to some Szego class measure, then asymptotic values of 
the monic recursion coefficients on any other finite support [a, h] are 



lim a' 



16 



Proof. It is evident that if © is satisfied, the linear transformation x — >■ mx + c 
which changes the measure support to [a, b] does not affect to its fulfilment since 
it is not singular in the interval, the rest is a simple application of the result 

m. □ 
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3 System-Reservoir structures 



The observables in an open quantum system are affected by the unavoidable 
interaction with the environment. This environment may be described by an 
infinite number (often called a reservoir) of bosonic or femionic modes labeled 
by some real number. The internal dynamics of the reservoir is given by some 
Hamiltonian of the form 

Hrcs = / dxg{x)alax, 



where in a physical context x could represent some continuous real variable such 
as the momentum of each mode, and Xmax tbe maximum value of it which is 
present in the reservoir (it could be infinity). In this picture g(x) represents 
the dispersion relation of the reservoir which relates the oscillator frequency 
to the variable x. The creation and annihilation operators satisfy the contin- 
uum bosonic [a^^jaj^] — S{x — y) or fermionic {ax,a^} — S{x — y) commutation 
rules. We assume that the frequencies g{x) and momenta x of the reservoir are 
bounded. 

The internal dynamics of the open quantum system are described by a un- 
specified local Hamiltonian operator H\oc and we assume that the interaction 
between the system and the reservoir is given by a linear coupling 

[ dxh{x)A{ax + al), 



where A is some operator of the system and the coupling (real) function h{x) 
describes the coupling strength with each mode. The total Hamiltonian for 
system plus reservoir is then given by 

H ^ Hioc + H^es + V = Hioc + / dxg{x)alax + / dxh{x)A{ax + al). 

Jo Jo 

(13) 

It has been shown that the dynamics induced in the quantum system by its 
interaction with the reservoir is completely determined by a positive function 
of the energy (or frequency uj) of the oscillators called the spectral density J{uj) 
[TJ[3]. For the continuum model of the reservoir we are considering, this function 
is given by, 

J(u;)=nh'[g'\u)]^^-^, (14) 
dcu 

where g~^[g{x)] — g[g~^{x)] = x. Physically, can be interpreted as 

the number of quanta with frequencies between lu and lo + Suj as Suj ^ i.e. it 
represents the density of states of the reservoir in frequency space. The spectral 
function thus describes the overall strength of the interaction of the system with 
the reservoir modes of frequency uj. 

This physical introduction motivate the next mathematical definition. 

Definition 3.1. We define the spectral density of a reservoir as a real non- 
negative and integrable function J{uj) inside of its (real positive) domain, which 
could be the entire half- line to G [0, oo). 
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Of course, given only a spectral density J{^jj), the dispersion relation g{x) 
and the coupling function h{x) are not uniquely defined. We shall make use of 
this freedom to implement a particularly simple transformation of the bosonic 
modes, and we will chose the dispersion function to be linear g{x) = gx. We 
now move on to our main theorem. 

Theorem 3.1. A system linearly coupled with a reservoir characterized by a 
spectral density J{uj) is unitarily equivalent to semi-infinite chain with only 
nearest-neighbors interactions, where the system only couples to the first site 
in the chain (see figure QP. In other words, there exist an unitary operator 
Un {x) such that the countably infinite set of new operators 

bl= f "^^^ dxUn{x)al (15) 



satisfy the corresponding commutation relations [6„,6tJ — 6nm for bosons, and 
{pmb]^} = 5nm for fcrmions, with transformed Hamiltonian 

H' = dx U„ {x)H = Hioc+CQA{bo + bl) + Y^ U„bibn+tnbi+ibn+tnbibn+l, 

(16) 

where cq, tn, w„ are real constants. 

Proof. The proof is by construction. Since J(w) is positive, h{x) is real, this 
defines the measure d^{x) = h^(x)dx. Then write 

U^{x) = h{x)p^{x) = h{x)^^, (17) 

where Pn{x) are some set of ortho normal polynomials with respect to the mea- 
sure dii,{x) = h^{x)dx with support on [0,Xmax]- Then it is clear that Un{x) is 
unitary (actually orthogonal as it is also a real transformation) in the sense of 



Jo 



dxUn{x)U*^{x) = / dxUn{x)Um{x) 

dxh^{x)pn{x)pm{x) = S„ 



SO the inverse transformation is just 

ai^J2Un{x)bl (18) 

n 

Moreover, for bosons 

[bn,bl^] = / dxdx'Un{x)Ujn{x')[ax,al,] 



dxdx Un{x)Um(x')5{x — x) 

dxUn{x)U„i{x) = Snm- 
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and similarly it is proved that {bn,bl^} — Snm for fermions. It remains to 
determine the structure of the transformed Hamiltonian H' , note that V is 
transformed like: 

V^AV] dxhix)U^ix)ib„ + hi) = i V Z""""" dxh^{x)^^{br. + bi), 

since for monic polynomials 7ro(a;) = 1 we find 

„ Jo hn\\ 

so Co = ||7ro||. For the H^cs term, note that with the choice of linear dispersion 
function g(x) = gx for the spectral density J(i^), one obtains 



= X! / dxg(x)Un{x)Urn{x)bl^bm 

= ^ / dxh'^{x)g{x)pn{x)pm{x)blbm 

= / dxxh^ {x)pn{x)pm(.x)b\j)m- 



With the recurrence relation @ we substitute the value of xpn{x) in the above 
integral to find 



„ .„ "'O 



1 A B ' 

—Pn+l{x) + y^Pnix) + -T^Pn-l{x) 



Pra{x)b\br, 



then orthonormality yields 

n 

where we have used the relation between monic and orthogomal recurrence 
coefficients So finally we have 



□ 

Remark 3.1. Note that from a topological point of view, mapping a continuous 
reservoir into a discrete chain is nothing but accounting for the separability of 
the topological space P (or the continuous functions space, if one prefers), and 
therefore the existence of a numerable basis. 
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Unix) 





Figure 1: Illustration of the effect of the transformation Un{x). On the left 
side system central is interacting with a reservoir with continuous bosonic or 
fermionic modes, paremetrized by x, after Un{x) the system is the first place of 
a discrete semi-infinite chain parametrized by n. 

Remark 3.2. In addition, the theorem provides us with a way to construct 
the exact mapping for any spectral density by simply looking for the family of 
orthogonal polynomials with respect to the measure d/x(x) = h?{x)dx. This can 
be done analytically for many important cases, as we demonstrate in section 
Ulwith the SBM, however even if the weight h'^{x) is a complicated function, 
families of orthogonal polynomials can be founded by using very stable numerical 
algorithms such as the ORTHPOL package |16) . 

Definition 3.2. We will say that a spectral density J(aj) belongs to the Szego 
class if the measure which it induced by d^ix) = h?{x)dx, with g{x) = gx, 
belongs to the Szego class. 

Corollary 3.1. For a Szego class spectral density J(cli) the tail of the semi- 
infinite chain tends to a translational invariant chain with 



Remark 3.3. Typical examples of Szego class spectral density are those strictly 
positive in its domain J{g{x)) > for x € [0, Xmax]- These includes a wide range 
of spectral densities with practical interest, particularly the ones used for the 
SBM discussed in the following section. 

Remark 3.4. The asymptotic properties of the frequencies and couplings for 
J{lj) in the Szego class have an important physical implication for open-quantum 
systems. Away from the quantum system, the chain becomes effectively trans- 
lationally invariant and it is a standard exercise in condensed matter physics 
to diagonalize such a chain in terms of plane waves |17| . For chains derived 
from a Szego class J{oj), the dispersion ujk of the eigenstates of the asymptotic 
region of the chain is uJk = Wc(l — cos(7rfc))/2 where the momentum runs from 
to 1 (cjc = 5(1))- This result shows rigorously that the asymptotic region 
can support excitations over the full spectral range of the original environment, 
as one would intuitively expect on physical grounds. The translational invari- 
ance of the asymptotic region means that excitations of the environment can 



lim Cli„ = 2 lim t„ = 



g-^mayi 



4 
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propagate away from the quantum system without scattering and are eflFectively 
"lost" irreversibly to the environment over time. One could also consider this 
asymptotic region of the chain as a secondary environment, in which case the 
system bath Hamiltonian can be recast as a finite chain in which the final mode 
is damped by the translationally invariant chain. In this picture the damping 
of the final member of the chain is universal for all spectral densities in the 
Szego class. Such a truncated picture might allow for improved efficiency in 
the numerical simulation of open-system dynamics, and the detailed physics be- 
hind this intuitively appealing picture of relaxation and dephasing in the chain 
picture will be dealt with elsewhere. 

Remark 3.5. A spectral function that appears in nature, most notably in 
photonic crystals, but which does not fall into the Szego class would be a spectral 
function that contains gaps, for instance 

J{ll)) — [Ji(lj)9{uJi — Uj) + J2{U})0{UJ — Ll!2)]0{uJc — 

Such a spectral function is zero in the region uji < uj < uj2. However, if Ji,2(w) 
are in the Szego class over [0, uji] and [w2, Wc] respectively, then the environment 
can be considered as two separate chains, for which the theorems above apply 
for each separate chain. 



4 Applications 

4.1 Continuous Spin-Boson spectral densities 

As a prototypical example let us consider the spin-boson Hamiltonian which 
describes the interaction of a TLS with an environment of harmonic oscillators 

mm 

HsB = Hioc + / dxg{x)alaj: + -0-3 / dxh{x){aj. + aj,), (19) 
where are bosonic operators, and the local Hamiltonian of the spin is given 

by, ^ ^ 

H\oc = + 2^^3- (20) 

where cti and ag, are the corresponding Pauli matrices. 

Firstly we will consider spectral functions bounded by a hard cut-off at an 
energy ujc, hence the cut-off Xmax = 5~^(wc) appearing in the integrals in ([T5|. 
which are usually parameterized as [S] 

J(a;) = 2Traujl-'Lj'0{uj - w^), (21) 

where a is the dimensionless coupling strength of the system-bath interaction 
and 6{uj — ujc) denotes the Heaviside step function. In the SBM literature, spec- 
tral functions with s > 1 are referred to as super-ohmic, s = 1 as Ohmic, and 
s < 1 as sub-Ohmic. The sub-ohmic and Ohmic cases are known to show inter- 
esting critical behaviour at zero temperaure as a function of a [H [3] • According 



13 



to (jl4p and our convention in taking linear dispersion relations, this continuous 
spectral function is related to the Hamiltonian parameters by, 



g{x) = LOcX, (22) 
h{x) = V2^LUcx'/'^. (23) 

With this choice of g{x) and h{x), Xmax = 1 and absorbing the common fac- 
tor uJc\^2a in the system operator A = the following matrix elements 
generate the mapping onto the chain, 

C/„(x)=a:^/2pi0'^)(x), (24) 

here P,^^'''\x) — Pn'^\x)N~^ is a (normalized) shifted Jacobi polynomial (from 
support [-1,1] to support [0,1]) for which the orthogonality condition reads 

dxx'P^^^'\x)P^^^'\x) = SnmK (25) 

K = } , (26) 



2n 



Therefore we get 



A j dxh{x){ax + aj.) = ^'^3 Y^(^o + ^1) 

where % = /g"" "/(w) du as defined by Bulla et al. in their NRG study of the 
SBM. They also use a mapping of the SBM onto a chain, however in their case 
the mapping is performed on an approximate discretised representation of the 
oscillator reservoir - see section 14.21 

For the transformed of Jq""^'' dxg{x)a\,ax we obtain from the recursion coef- 
ficients of the Jacobi polynomials !^18^ , frequencies and tunneling amplitudes to 
be 

^ Y(^+(. + 2n)(2 + s + 2n))' 



Wc(l + n)(l + s-f n) l^ + s + 2n ^^g) 



{s + 2 + 27i)(3 + s -f 2n) V l + s + 2n 

It is quite easy to check that the measure dfi{x) = x'^dx belongs to the Szego 
class for any s < oo, the integral (O on [0, 1] is 

\ii{x'^)dx , 
^ ^ - -7rsln(2). 



^l-{2x-lY 

1 ^1 

Thus lim„^oo = 5 ^^'^ lim„->oo /^n = ^ = T6' 

hm iLj„ = — , hm t„ = — 

n— >-oo Z n— >oo 4 



as it can be seen directly from the expressions (j27p and ((28|) . Interestingly, 
smaller values of s lead to faster convergence of a;„ , i„ towards these limits as n 
increases. 
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If instead of a hard cut-off we use an exponential cut-off for the spectral 
function 

J{lu) ~ 2TiaLo\~^uj^e~'^ , 
and an integral on [0, oo), then for g{x) — ujcX we find 

h{x) = V2^u:cx''^e-'/^, 
the transformation to the chain can achieved using 

[/„(x) = x^'^UX^) = x^'^e-^'^Ll,{x)N~\ (29) 

where L!^{x) is an associated Laguerre polynomial and N„ is the normalisation 
of the orthogonality relation for the associated Laguerre polynomials: 

1 

dxx'e~''L'^{x)Ll{x) = (5„™iV2 (gg^ 

N?. = l^^^±I±^. (31) 
n! 

So the transformation reads 

/■°° 1 
A J dxh{x){ax + al) = -0-300(60 + 6^), 







with Co = LjJc\j2aT{s + 1) = where 770 = Jj, J{uj)duj\ and for the other 
term from the recurrence relations of the associated Laguerre polynomials |18) , 
the frequencies and tunneling are 

= (2n + 1 + s) , (32) 



tn = ujc\/{n + l){n + s + l).. (33) 

4.2 Discrete spectral densities - The logarithmically dis- 
cretized Spin-Boson model 

Our main result (Theorem I3.ip maps continuous environments onto discrete 
semi-infinite chains which are are ideal for t-DMRG simulation. However in 
some situations we may want to simulate systems which are coupled to a discrete 
set of oscillators which are not described by a continuous spectral function. This 
may occur physically in many ways, for example systems coupled to vibrations 
of small or finite-sized environments are coupled to discrete spectral densities. 
Many of these discrete environments can be still be transformed into a chain 
using the same procedure we have presented for the continuous mapping, the 
only different is that one now makes use of the classical discrete orthogonal 
polynomials [TO]. One example of this is the finite-size discretisation of the 
spectral densities considered in section 14.11 It will be shown in section 14.31 
that these discrete environments can be mapped onto chains suitable for t- 
DMRG analysis using the classical Hahn polynomials. In this subsection we 
wish to focus on another important class of discrete spectral densities which are 
logarithimcally discretised approximations to continuous environments. These 
discretised spectral densities are a key part of the powerful NRG algorithm used 
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to study condensed matter systems, and a detailed description of this numerical 
method and the reasons for performing this discretisation are given in Ref. [S]- 
In the NRG approach to the same SBM considered in section 21 the con- 
tinuous spin-boson Hamiltonian is approximated by representing all the modes 
in a energy range WcA"" — WcA"""^ by a single effective degree of freedom 
described by creation and annihilation operators a„,ajj. These modes may 
either be fermionic or bosonic. The logarithmic discretisation parameter A 
is always greater than unity, and the original continuum spin-boson Hamilto- 
nian is recovered by sending A — 1. Assuming a spectral density of the form 
J(a;) = 2Tr au!l~^uj^ 6 (lUc—lo), this logarithmic discretisation procedure generates 
a discrete spin-boson Hamiltonian Hl of the form [9], 

oo oo 

Hl = Hioc + ^ 7n(a« + ai) + ^ Cnaia„, (34) 

V n=0 n=0 



where 



^2 ^ ^^2(l_^-{i+s))^-n(l+s) ^^^^-Hl+s)^ ^3J3^ 
1 + S 

Cn = — ^ — — w,A-" = GA-". (36) 

As described in detail in Ref. W , the NRG algorithm can only be applied after 
the Hamiltonian Hl of Eq. (IMl) is mapped onto a. 1 — D nearest-neighbour 
chain of the form, 

+ t,,bibn+l. (37) 



To do this Bulla et al. [5] first introduce a new set of bosonic modes 6- 
defined by an orthogonal transformation U which acts on the a„,ajj, 

bi = (38) 

m 

bn — ^ ^ t-^nm^m, ■ (^'^) 
m 

They then derive the following recurrence relation for the matrix elements Unm 
that bring the Hamiltonian of to the chain form of p7p [5], 

CnUmn ^mUmn + ^mUm+ln + ^m — l^-^m— In- (^0) 

This recurrence relation is then solved numerically by an iterative procedure 
to provide the hopping parameters i„ and site energies ojn which are the input 
parameters for the NRG algorithm. This recurrence relation is the standard 
method for finding the chain parameters for almost all numerical simulation 
methods that use the chain representation for impurity problems at present. 
However the numerical solutions of the recurrence relations are highly sensitive 
to numerical noise and are often unstable as the dimension of U becomes large. 
In Ref. fS] it is stated that numerical instabilities typically appear after about 
25 — 30 iterations. The results presented in this paper show that the problem 
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of implementing this chain mapping effectively amount to the problem of find- 
ing a suitable basis of orthogonal polynomials to represent the bath degrees of 
freedom. The recursive method given by Bulla et al. [5] is in fact equivalent 
to the recursive Stieltjes method for determining the coefficients of orthogonal 
polynomials, a procedure which is well-known to be numerically unstable for 
many weight functions, and in certain cases is so inaccurate that it completely 
fails for all practical purposes [30]. As orthogonal polynomials are very useful 
and widely used in numerical quadrature, the problem of determining recurrence 
coefficients accurately for arbitrary weight functions has been investigated in de- 
tail, and alternative numerical algorithms for determining recurrence coefficients 
to arbitrary accuracy have been developed, most notably the ORTHPOL pack- 
age of Gautschi [IB]. Our results showing that the chain mapping can always 
be cast as a problem of finding recurrence coefficients of orthogonal polynomials 
has thus allowed us to make use of these methods, therefore removing the issue 
of numerical stability from the general problem. 

However, for the specific example considered in this section we can in fact 
solve the recursion relation analytically, and can give explicit formulae for the 
parameters i„ , aj„ . We now demonstrate that the solution Unm to the recurrence 
relation of [9 can be written as, 

_ A~^P„(A-",A-M|A-i) 

where the functions p„(A"'", A"^ l|A"i) in Eq. gl]) are known as the little-q 
Jacobi polynomials [21]. These polynomials obey the following orthogonality 
relation which defines the normalization constants iV„ that appear in Eq. (|4ip , 

oo 

SnmN^ = ^A-'=(l+^)p„(A-^A-^l|A-l)p„(A-^A-^l|A-l) (42) 

A;=0 

^ - (43) 

(A-(''+i); A-i)2(l - A-(2"+i+''))' ^ ' 

In Eq. (a : q)n is the q-shifted factorial described in [3T], and defined as 

(a : q)n = (1 - a)(l - aq){l - aq^)...{l - aq"-^). (44) 
Noting that Pq{x, a, b\q) = 1 [3T], 

oo 

^A-^■(l+V(A-^A-^l|A-l) = 6onN^ (45) 

fe=0 

= ^o» ^_^i(,+i) - (46) 

One further property of the little-q Jacobi polynomials that we require to per- 
form the chain mapping is their three-term recurrence relation. This is given 

by, 

A->,(A-",A-M|A-i) = (A,-|-Q)p,(A-",A-M|A-i) 

- A,p,+i(A-",A-M|A-i) 

- Qp,_l(A-^A-^l|A-l), (47) 
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where, 

(1 - A-(2J + 1+''))(1 - A-(2j+2+s))' 

(1- A-J)2 
(1 - A-(2J+''))(1 - A-(2j+i+'s 

Using the orthogonahty relations one can easily show that the transformation 
U is an orthogonal matrix and that the 6„ , 6jj operators obey the same bosonic 
or fermionic commutation relations as the a„,ajj operators. 

Now we express the original Hamiltonian in terms of the 6n,^l operators. 
The second term of (l34ll transforms as, 



oo oo 



n=0 n=Ok=Q 



= 7. E ^k'ih + bl) A-"(i+^V(A-", A-^ l|A-i) 

oo 
k=Q 



= jsNoibo + bl) = ^]{bn + bl) 



where we have used property of Eq. (|46p and in the last step the explicit form 
of JsNq has been used, 

where 770 is the coupling defined by Bulla et al. [3]. The chain couplings and 
energies arise from the transformation of the third term of which then 

becomes, 

00 00 00 



00 



X A-"(i+'*) A->j(A-", A-^ l|A-i) pfe(A-", A-^ l|A-i). 

n=0 

The orthogonality condition cannot be used directly here as the weight function 
of the sum A^"^'*^^) is multiplied by an extra factor of A^". However, we 
can absorb this factor in the sum using the recurrence relation of the little-q 
polynomials on the braced term above. Substituting the RHS of the recurrence 
relation into Eq. ([501 allows the orthogonality relations to be used, and ensures 
that the resulting couplings are nearest-neighbour only. This generates a 1 — D 
chain Hamiltonian H' of the form of (1161) . with energies w„ and hoppings i„ 
given by, 

w„ = Qs{An + Cn), (50) 

tn = -Cs(^.)a„ (51) 
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We now show that our exphcit form for Umn in terms of httle-q Jacobi 
polynomials satisfies the recurrence equation of Bulla et al. Eq. PO)) . We first 
use Eqs. (15(11) and ([51]) to express the RHS of (^0]) in terms of the little-q Jacobi 
recurrence coefficients yl„,C„. Substituting C„ — CsA^" on the LHS as well, 
one obtains 

~ Cs ( ) A,n-lUm-ln- (53) 

Now substitute for the Umn using (j41l) we obtain, 

A->™(A-^A-^l|A-l) - (A„ + C™)p„(A-",A-M|A-i) 



A„-lP™-l(A-",A-^l|A-l) 



+ AnP„+l(A-'^A-^l|A-l). 

The final additional that result we need to use is that the following relations 
holds. 

Cm - Am-1 f ™ > 1, (54) 

Co = m = 0, (55) 

which is easily shown from the definitions of Am^Cm- Substituting this into 
Ea. ((56l) . the equation reduces to, 

A->„(A-",A-M|A-i) - (A„ + C„)p™(A-",A-M|A-i) 

+ A„p™+l(A-^A-^l|A-l) 
+ C„p„,_l(A-^A-^l|A-l). 

This is just the recurrence relation for the little-q Jacobi polynomials. Therefore 
the original Unm of Eq. (HT|) is the solution of the recurrence relations of Bulla 
et al.. The energies and hopping amplitudes are simple functions of An, Cn and 
the following additional relation can also be derived using Eqs. (|5n|) and ([5T|) . 



Nn+lJ \ Nn 

4.3 Linearly-discretised baths 

We now consider an environment consisting in a large, but discrete number 
+ 1 of environmental degrees of freedom. The Hamiltonian for such a system 
can be written as 

N N 

H ^ H10C + YJ^^ln{an + al) + ^C,nalan, (56) 

^''^ n=0 n=0 

where C„ = LUcn/{N + 1), the coupling constants 7„ are 
2 o f s + n 



In = , (57) 
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and the bosonic operators obey the discrete commutation relation [a„,aj„] = 
Snm- The expressions for C„ and 7„ arise from a hnear discretisation of the 
continuous spectral function J{(jj) = 2naujl~''uj'' with N +1 oscillators equally- 
spaced in k approximating the bath. As TV — ?> cxd we recover the continuum limit 
which was mapped using Jacobi polynomials. This discrete Hamiltonian can be 
mapped to the nearest-neighbour chain using the orthogonal transformation 



N 

bk = ^i ] Qk{n,s)p-^an, (58) 



n=0 



where Qi{n,s) — Qi{n,s,0,N) is a Hahn polynomial [HI [51]. The Hahn poly- 
nomials are classical discrete orthogonal polynomials and are defined by the 
discrete inner product 

E ( "t"" ) ( ^%^_l'')Qiin,a,f3,N)Q,{n,a,f3,N)^plS,,. (59) 

Introducing the Pochhammer symbol {a)k — a{a + l)...{a + N — 1), the normal- 
ization factors pk are given by 

^ (-l)^(fc + a + /3 + l)jv+i(/3 + l)fcfc! 

(2fc + a + /3 + l)(a + l)fe(-7V)feiV! • ^"""^ 

Using the inner product relation is it easy to show that the new modes bk, bj, obey 
the same commutation relations as the original modes, and using the identity 

s + x\/' 5 + y \ Qk{x,a,l3,N)Qkiy,a,P,N) ^ ^ 
^ / V y / ^0 Pk -v' ^ ) 

one can show that the inverse transform is 

N 



an = 

k=0 



Y.[ Z ] Q'^i^^^)h- (62) 



The Hahn polynomials also obey the three-term recurrence relation 

xQnix,a,l3,N) = A„Qn+iix,a,l3,N) - {An + Cn)Qn{x,a,l3,N) 

+ CnQn-i{x,a,(3,N), (63) 



where 



^ {n + a + /3 + l)(n + a + 1){N - n) 

i2n + a + /3){2n + a + f3 + l) ' ^ ' 

^ ^ n{n + a + P + N+l){n + P) 

{2n + a + /3){2n + a + /3 + 1)' ^ ' 



We now substitute (1621) into H to obtain the chain Hamiltonian. As in the 
previous examples, the three-term recurrence and orthogonality of the Qk(n, s) 
ensures the generation of a discrete nearest-neighour chain with energies and 
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couplings determined by the recurrence coefficients of tlie Hafin polynomials. 
The transformed Hamiltonian He is given by, 



c 



Hloc + 7: 




(7z{ha + fej) + X! + tnhlj)n+l + tn^l+i&n, (66) 



n=0 



'n 




(68) 



"tl 



N + l p, 



The discrete weight function is not in the Szego class, and at finite N we 
find that w„ — >■ uJc/4:,tn as n ^ N. However, taking the limit N ^ oo 
leads to the recovery of the results obtained for a continuous power-law spectral 
density in section HTTl Formally this results from the limit relation |21) . 



although this can also be seen by using taking the limit of the explicit expres- 
sions for the recurrence relations An , C„ and normalizations Nn for the Hahn 
polynomials. Using the explicit expressions, it is also easy to show that the dif- 
ference between the chain energies and hoppings for the continuous and linearly 
discretized bath are of 0{n/N) for n <^ N. The differences in the behaviour of 
the chain parameters in the finite and continuous cases has a physical signifi- 
cance that will be explored in detail elsewhere. Here, it will be enough to say 
the vanishing of the t„ as n — > oo is related to the existence of a recurrence time 
in the dynamics of the discrete problem. 

5 Conclusions 

In this article we have demonstrated how the problem of a quantum system cou- 
pled linearly to the coordinates of a continuous bosonic or fermionic environment 
can be mapped exactly onto a one-dimensional model with only short-range in- 
teractions. Using the theory of orthogonal polynomials, we have shown how 
the parameters of this one dimensional representation can be simply obtained 
from the recurrence coefficients of the orthogonal polynomials of the spectral 
function describing the environment. Certain important spectral functions cor- 
respond to the weight functions of the classical orthogonal polynomials, and 
the parameters of the chain Hamiltonian can be given in closed analytical form. 
When the spectral function is not a classical weight function, we can make use 
of numerical techniques developed for determining recurrence coefficients that 
allow us to find the chain parameters to high precision and for arbitrary long 
chains. Our use of orthogonal polynomials simultaneously removes the need to 
discretise the environmental degrees of freedom and also the numerical instabil- 
ities that plague the recursive method previously used to implement this chain 
transformation. This is particularly important, as hitherto the principal use of 
the chain transformation is to enable numerical simulation of the dynamics and 
ground state properties of systems which strongly interact with their environ- 
ments. The work presented here was principally motivated by an attempt to 
use the powerful t-DMRG simulation technique to study strongly-interacting 



hm QniNx,a,b,N) = {-irP,-/i2x-l)/P.^,^\l), 
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open-quantum systems, a problem of considerable topically, especially in the 
emerging field of quantum effects in biology. The combination of our exact 
mapping has now made it possible to use t-DMRG to obtain extremely accu- 
rate simulations of open-system dynamics with no artefacts arising from having 
a discretised representation of the environment and a numerical approximation 
to the mapping. Using the theory of orthogonal polynomials, we have also given 
examples that could be used in the NRG community. 

In addition to these practical advantages, the analytical mapping reveals an 
interesting structure to the system-environment systems we consider, in particu- 
lar, the result that the chain representations of all spectral functions in the phys- 
ically broad, Szego class look the same in the asymptotic limit. This observation 
may lead to further improvements in the numerical efficiency of open-quantum 
simulations, allowing the chain representations to be truncated in an efficient 
way that is universal for all spectral functions in the Szego class. We have 
also presented results showing the relations between the chain representations 
of linearly and logarithmically discretised spin-boson models and the continuous 
case, results which may also prove useful in determining how to optimize the 
accuracy of simulations with a discretised representations. 
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